Data inladen

Baits<-read.table(file = "~/GitHub/vespa_analyses/Input/Baits_useful.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))

Flights<-read.table(file = "~/GitHub/vespa_analyses/Input/Flights_KMI.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))

Nests<-read.table(file = "~/GitHub/vespa_analyses/Input/Nests_useful.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))

Shortind_incomplete<-read.table(file = "~/GitHub/vespa_analyses/Input/Shortind_incomplete.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))

Shortpot_incomplete<-read.table(file = "~/GitHub/vespa_analyses/Input/Shortpot_incomplete.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))

Individuals<-read.table(file = "~/GitHub/vespa_analyses/Input/Individuals.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))

Aanpassingen datasets

# Overbodige kolommen wegdoen, foutje in excel
Flights<-Flights[1:(length(Flights)-1)]
Shortind_incomplete<-Shortind_incomplete[1:(length(Shortind_incomplete)-2)]
Shortpot_incomplete<-Shortpot_incomplete[1:(length(Shortpot_incomplete)-2)]

#Foerageersnelheid en vliegsnelheid (m/s) toevoegen aan tabel
Flights$test_rule_of_thumb<-Flights$Distance/(Flights$Flighttime_min*60)
Shortind_incomplete$test_rule_of_thumb<-Shortind_incomplete$Distance/(Shortind_incomplete$Flighttime_min*60)
Shortpot_incomplete$test_rule_of_thumb<-Shortpot_incomplete$Distance/(Shortpot_incomplete$Flighttime_min*60)

Flights$ForagingSpeed<-Flights$Distance*2/(Flights$Flighttime_min*60)
Shortind_incomplete$ForagingSpeed<-Shortind_incomplete$Distance*2/(Shortind_incomplete$Flighttime_min*60)
Shortpot_incomplete$ForagingSpeed<-Shortpot_incomplete$Distance*2/(Shortpot_incomplete$Flighttime_min*60)


## Urbanisatiecirkels per vliegtijd toevoegen
Flights<-merge(Flights, Baits[, c("NestID", "BaitID", "Urbanisation25m", "Urbanisation50m", "Urbanisation100m")], by=c("NestID", "BaitID"), all=TRUE)
Shortind_incomplete<-merge(Shortind_incomplete, Baits[, c("NestID", "BaitID", "Urbanisation25m", "Urbanisation50m", "Urbanisation100m")], by=c("NestID", "BaitID"))
Shortpot_incomplete<-merge(Shortpot_incomplete, Baits[, c("NestID", "BaitID", "Urbanisation25m", "Urbanisation50m", "Urbanisation100m")], by=c("NestID", "BaitID"))

## Urbanisatietrajecten per vliegtijd toevoegen
Flights<-merge(Flights, Baits[, c("NestID", "BaitID", "Traject25m", "Traject50m", "Traject100m")], by=c("NestID", "BaitID"), all=TRUE)
Shortind_incomplete<-merge(Shortind_incomplete, Baits[, c("NestID", "BaitID", "Traject25m", "Traject50m", "Traject100m")], by=c("NestID", "BaitID"))
Shortpot_incomplete<-merge(Shortpot_incomplete, Baits[, c("NestID", "BaitID", "Traject25m", "Traject50m", "Traject100m")], by=c("NestID", "BaitID"))

## Gewicht individu per vliegtijd toevoegen
Flights<-merge(Flights, Individuals[, c("NestID", "BaitID","ColorInd", "Weight_ind")], by=c("NestID", "BaitID", "ColorInd"), all=TRUE)
Shortind_incomplete<-merge(Shortind_incomplete, Individuals[, c("NestID", "BaitID", "ColorInd", "Weight_ind")], by=c("NestID", "BaitID", "ColorInd"))
Shortpot_incomplete<-merge(Shortpot_incomplete, Individuals[, c("NestID", "BaitID", "ColorInd","Weight_ind")], by=c("NestID", "BaitID", "ColorInd"))

## Nesthoogte per vliegtijd toevoegen
Flights<-merge(Flights, Nests[, c("NestID", "Height")], by=c("NestID"), all=TRUE)
Shortind_incomplete<-merge(Shortind_incomplete, Nests[, c("NestID", "Height")], by=c("NestID"))
Shortpot_incomplete<-merge(Shortpot_incomplete, Nests[, c("NestID", "Height")], by=c("NestID"))

library(ggplot2)
library(ggpubr)

Volgende datasets zijn gecreëerd:

Voor het model: Flight time ~ Distance gebruiken we de dataset vliegtijden per individu Omdat we hiermee de theoretische regel 1min=100m kunnen verifiëren. De imkers nemen hiervoor altijd kortste meting

Voor modellen met weerparameters, gewicht, urbanisatie: Hiervoor nemen we telkens de hele dataset, omdat elke meting van deze factoren afhangt.

Datapunten op 2km weglaten

We gaan verder met de dataset zonder de outliers in Melle

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
Flights2<-Flights %>%  filter(NestID!=1)
Shortind2_incomplete<-Shortind_incomplete %>%  filter(NestID!=1)
Shortpot2_incomplete<-Shortpot_incomplete %>%  filter(NestID!=1)

# ID maken per individu
Flights2$Individualcode <- paste(Flights2$NestID, Flights2$BaitID,Flights2$ColorInd, sep="_")
Shortind2_incomplete$Individualcode <- paste(Shortind2_incomplete$NestID, Shortind2_incomplete$BaitID,Shortind2_incomplete$ColorInd, sep="_")

VRAAG: Hoe grafiek maken van parameters mixed model?

1.1) ForagingSpeed vs Temperature

Graphs

Model Output

  • Normality niet ok!

  • Niet significant!

library(lmerTest)
## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(lme4)

model_temp_all<-lmer(ForagingSpeed ~ Temperature_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model_temp_all)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## ForagingSpeed ~ Temperature_KMI + (1 | NestID) + (1 | NestID:BaitID) +  
##     (1 | Individualcode)
##    Data: Flights2
## 
## REML criterion at convergence: 595.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.19526 -0.45657  0.04751  0.38435  3.07187 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1802   0.4245  
##  NestID:BaitID  (Intercept) 0.7973   0.8929  
##  NestID         (Intercept) 0.7893   0.8884  
##  Residual                   0.3540   0.5949  
## Number of obs: 230, groups:  Individualcode, 91; NestID:BaitID, 68; NestID, 34
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       3.234075   0.537550 135.595552   6.016 1.57e-08 ***
## Temperature_KMI  -0.001193   0.026258 175.559355  -0.045    0.964    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temprtr_KMI -0.919
anova(model_temp_all, ddf="Satterthwaite", type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq    Mean Sq NumDF  DenDF F value Pr(>F)
## Temperature_KMI 0.00073042 0.00073042     1 175.56  0.0021 0.9638
res<-residuals(model_temp_all)
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.98221, p-value = 0.005473
qqnorm(res)
qqline(res)

1.2) Flight time vs Distance + Temperature

Graph

Klopt dit? Kan dit ook met mixed model?

library(knitr)
library(kableExtra)
## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 2 > 1' in coercion to 'logical(1)'
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library("plot3D")

x <- Flights2$Distance
y <- Flights2$Temperature_KMI
z <- Flights2$Flighttime_min

fit <- lm(z ~ x + y, na.action=na.exclude)
x.pred <- seq(min(x[!is.na(x)]), max(x[!is.na(x)]), length.out = 20)
y.pred <- seq(min(y[!is.na(y)]), max(y[!is.na(y)]), length.out = 20)
xy <- expand.grid( x = x.pred, y = y.pred)
z.pred <- matrix(predict(fit, newdata = xy), 
                 nrow = 20, ncol = 20)
fitpoints <- predict(fit)

scatter3D(x, y, z, pch = 19, cex = 0.6, colvar=FALSE, col="dodgerblue3", theta = 210, phi = 10, bty="u", col.panel ="grey93", expand =0.4, col.grid = "white", xlab = "Distance", ylab = "Temperature", zlab = "Flight time", surf = list(x = x.pred, y = y.pred, z = z.pred,  
facets = TRUE, col=ramp.col(col = c("dodgerblue4", "seagreen2"), n = 100, alpha=0.8), fit = fitpoints, border="black"),main = "Flight time vs Distance + Temperature")

Met plotly

Had met deze link problemen met expand.grid https://stackoverflow.com/questions/38331198/add-regression-plane-to-3d-scatter-plot-in-plotly Dan maar grid dat van hierboven gebruikt, maar nu klopt er precies iets niet? Alle datapunten liggen boven het oppervlak.

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(tidyverse)
## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──
## ✔ tibble  3.1.8     ✔ purrr   1.0.1
## ✔ tidyr   1.2.1     ✔ stringr 1.5.0
## ✔ readr   2.1.4     ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand()          masks Matrix::expand()
## ✖ plotly::filter()         masks dplyr::filter(), stats::filter()
## ✖ kableExtra::group_rows() masks dplyr::group_rows()
## ✖ dplyr::lag()             masks stats::lag()
## ✖ tidyr::pack()            masks Matrix::pack()
## ✖ tidyr::unpack()          masks Matrix::unpack()
FlightsNoNA<-Flights2 %>%  filter(Flighttime_min!="NA")
FlightsNoNA<-FlightsNoNA %>%  filter(Temperature_KMI!="NA")

fig <- plot_ly(FlightsNoNA, x = ~Distance, y = ~Temperature_KMI, z = ~Flighttime_min, size=1)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = 'Distance'),
                     yaxis = list(title = 'Temperature_KMI'),
                     zaxis = list(title = 'Flight time (min)')))
p2 <- add_trace(p = fig,
                z = z.pred,
                x = seq(2100, 0, by = -100),
                y = seq(0, 30, by = 10),
                type = "surface")
p2

Model Ouput

  • Normality niet ok!

  • Niet significant

model_tempdist_all<-lmer(Flighttime_min ~ Distance + Temperature_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.exclude, data=Flights2)
## boundary (singular) fit: see help('isSingular')
summary(model_tempdist_all)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Distance + Temperature_KMI + (1 | NestID) +  
##     (1 | NestID:BaitID) + (1 | Individualcode)
##    Data: Flights2
## 
## REML criterion at convergence: 1005.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2425 -0.4262 -0.1822  0.1635  4.4765 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1443   0.3799  
##  NestID:BaitID  (Intercept) 1.8910   1.3751  
##  NestID         (Intercept) 0.0000   0.0000  
##  Residual                   3.2814   1.8115  
## Number of obs: 230, groups:  Individualcode, 91; NestID:BaitID, 68; NestID, 34
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      3.5399788  0.9074782 66.9650223   3.901 0.000225 ***
## Distance         0.0056150  0.0008062 62.3323406   6.965 2.42e-09 ***
## Temperature_KMI -0.0686672  0.0445457 91.2883267  -1.541 0.126655    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Distnc
## Distance    -0.336       
## Temprtr_KMI -0.903 -0.010
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
anova(model_tempdist_all, ddf="Satterthwaite", type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Distance        159.173 159.173     1 62.332 48.5081 2.415e-09 ***
## Temperature_KMI   7.797   7.797     1 91.288  2.3762    0.1267    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res<-residuals(model_tempdist_all)
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.7916, p-value < 2.2e-16
qqnorm(res)
qqline(res)

2) ForagingSpeed vs Cloudcoverage

Graphs

ggplot(Flights2, aes(x=Cloudcoverage_KMI, y=ForagingSpeed)) + geom_point(color="lightblue") + geom_smooth(method="lm", formula =y ~ x, se=TRUE, fullrange=FALSE, level=0.95, color="darkblue") + ggtitle("All data KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))

Model Ouput

  • Licht significant

  • Normality niet ok!

model_cloud_all<-lmer(ForagingSpeed ~ Cloudcoverage_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model_cloud_all)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## ForagingSpeed ~ Cloudcoverage_KMI + (1 | NestID) + (1 | NestID:BaitID) +  
##     (1 | Individualcode)
##    Data: Flights2
## 
## REML criterion at convergence: 491.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2243 -0.5101  0.0644  0.4260  2.9113 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1866   0.4320  
##  NestID:BaitID  (Intercept) 0.9826   0.9912  
##  NestID         (Intercept) 0.7492   0.8655  
##  Residual                   0.3547   0.5956  
## Number of obs: 191, groups:  Individualcode, 74; NestID:BaitID, 52; NestID, 23
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)         3.2816     0.2715  24.9103  12.088 6.45e-12 ***
## Cloudcoverage_KMI  -0.4237     0.2090 141.7036  -2.027   0.0445 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Cldcvrg_KMI -0.326
anova(model_cloud_all, ddf="Satterthwaite", type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Cloudcoverage_KMI 1.4581  1.4581     1 141.7  4.1106 0.04449 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res<-residuals(model_cloud_all)
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.98299, p-value = 0.0204
qqnorm(res)
qqline(res)

3.1) ForagingSpeed vs Windspeed

Graphs

Telkens voor de modellen:

ForagingSpeed ~ Windspeed

ForagingSpeed ~ Windspeed²

plot19<-ggplot(Flights2, aes(x=WindSpeed_KMI, y=ForagingSpeed)) + geom_point(color="khaki3") + geom_smooth(method="lm", formula =y ~ x, se=TRUE, fullrange=FALSE, level=0.95, color="olivedrab") + ggtitle("ForagingSpeed ~ Windspeed  KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))

plot22<-ggplot(Flights2, aes(x=WindSpeed_KMI, y=ForagingSpeed)) + geom_point(color="khaki3") + geom_smooth(method="lm", formula =y ~ I(x^2), se=TRUE, fullrange=FALSE, level=0.95, color="olivedrab") + ggtitle("ForagingSpeed ~ Windspeed²  KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))


ggarrange(plot19, plot22 + rremove("x.text"), 
          labels = c("A", "B"),
          ncol = 2, nrow = 1)

Model Output

  • Beide modellen hoogsignificant

  • Normality voor beiden ok!

model_wind_all<-lmer(ForagingSpeed ~ WindSpeed_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model_wind_all)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ WindSpeed_KMI + (1 | NestID) + (1 | NestID:BaitID) +  
##     (1 | Individualcode)
##    Data: Flights2
## 
## REML criterion at convergence: 572.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.55366 -0.52807  0.07123  0.41517  2.99472 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1681   0.4100  
##  NestID:BaitID  (Intercept) 0.7091   0.8421  
##  NestID         (Intercept) 0.8242   0.9079  
##  Residual                   0.3160   0.5621  
## Number of obs: 230, groups:  Individualcode, 91; NestID:BaitID, 68; NestID, 34
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     3.94765    0.25805  70.73569   15.30  < 2e-16 ***
## WindSpeed_KMI  -0.19399    0.03984 165.66865   -4.87 2.59e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.584
anova(model_wind_all, ddf="Satterthwaite", type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## WindSpeed_KMI 7.4929  7.4929     1 165.67  23.713 2.594e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model_wind2_all<-lmer(ForagingSpeed ~ I(WindSpeed_KMI^2) + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model_wind2_all)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## ForagingSpeed ~ I(WindSpeed_KMI^2) + (1 | NestID) + (1 | NestID:BaitID) +  
##     (1 | Individualcode)
##    Data: Flights2
## 
## REML criterion at convergence: 580
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.26566 -0.54124  0.08453  0.41515  2.93663 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1414   0.3760  
##  NestID:BaitID  (Intercept) 0.7246   0.8513  
##  NestID         (Intercept) 0.8080   0.8989  
##  Residual                   0.3295   0.5740  
## Number of obs: 230, groups:  Individualcode, 91; NestID:BaitID, 68; NestID, 34
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          3.599231   0.226024  45.995965   15.92  < 2e-16 ***
## I(WindSpeed_KMI^2)  -0.022824   0.005175 191.376210   -4.41 1.72e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## I(WS_KMI^2) -0.387
anova(model_wind2_all, ddf="Satterthwaite", type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## I(WindSpeed_KMI^2) 6.4084  6.4084     1 191.38  19.451 1.719e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res<-residuals(model_wind_all)
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.98923, p-value = 0.0834
qqnorm(res)
qqline(res)

res<-residuals(model_wind2_all)
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.98929, p-value = 0.08548
qqnorm(res)
qqline(res)

3.2) ForagingSpeed vs Windspeed | afhankelijk van de windrichting

Hiervoor werd telkens voor elke meting nagegaan of de hoornaar met meewind (tailwind), tegenwind (upwind) of loodrechte wind (perpendicular) te maken had. Dit volgens de formules:

|𝜃flight −𝜃wind| ≤ 45 is tailwind

45 < |𝜃flight −𝜃wind|<135 is (quasi) perpendicular

|𝜃 flight −𝜃wind| ≥135 upwind

Graphs

ggplot(data=subset(Flights2, !is.na(Wind)), aes(x= Wind, col=Wind, y=ForagingSpeed)) + geom_boxplot()
## Warning: Removed 14 rows containing non-finite values (`stat_boxplot()`).

ggplot(data=subset(Flights2, !is.na(Wind)), aes(x= WindSpeed_KMI, col=Wind, y=ForagingSpeed)) + geom_point() + facet_grid(~Wind) + geom_smooth(method="lm", formula = y ~ I(x^2), se=TRUE, fullrange=FALSE, level=0.95)
## Warning: Removed 14 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 14 rows containing missing values (`geom_point()`).

Model Output

Anova van Foragingspeed en de 3 windinvloeden

  • interpretatie?
windanova<-lmer(ForagingSpeed ~ Wind + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.exclude, data=Flights2)
summary(windanova)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ Wind + (1 | NestID) + (1 | NestID:BaitID) + (1 |  
##     Individualcode)
##    Data: Flights2
## 
## REML criterion at convergence: 592.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1864 -0.4629  0.0327  0.3908  3.0729 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1877   0.4332  
##  NestID:BaitID  (Intercept) 0.7814   0.8840  
##  NestID         (Intercept) 0.7900   0.8888  
##  Residual                   0.3548   0.5956  
## Number of obs: 230, groups:  Individualcode, 91; NestID:BaitID, 68; NestID, 34
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    3.24826    0.22447  40.83023  14.471   <2e-16 ***
## Windtailwind  -0.14580    0.19936 214.92685  -0.731    0.465    
## Windupwind    -0.02152    0.18402 182.45501  -0.117    0.907    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Wndtlw
## Windtailwnd -0.240       
## Windupwind  -0.279  0.224

ForagingSpeed vs WindSpeed voor Tailwind dataset

  • Hoogsignificant

  • Normality ok!

data_tailwind<-Flights2 %>%  filter(Wind!="tailwind")

model_tailwind<-lmer(ForagingSpeed ~ WindSpeed_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=data_tailwind)
summary(model_tailwind)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ WindSpeed_KMI + (1 | NestID) + (1 | NestID:BaitID) +  
##     (1 | Individualcode)
##    Data: data_tailwind
## 
## REML criterion at convergence: 476.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.45714 -0.55810  0.04615  0.41914  2.91009 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1738   0.4169  
##  NestID:BaitID  (Intercept) 0.7757   0.8807  
##  NestID         (Intercept) 0.7954   0.8919  
##  Residual                   0.3318   0.5760  
## Number of obs: 188, groups:  Individualcode, 75; NestID:BaitID, 54; NestID, 31
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     3.81145    0.27541  60.89085  13.839  < 2e-16 ***
## WindSpeed_KMI  -0.16562    0.04216 130.28317  -3.929 0.000138 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.574
anova(model_tailwind, ddf="Satterthwaite", type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## WindSpeed_KMI  5.121   5.121     1 130.28  15.434 0.0001379 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res<-residuals(model_tailwind)
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.98883, p-value = 0.1474
qqnorm(res)
qqline(res)

ForagingSpeed vs WindSpeed voor Perpendicular dataset

  • Niet significant

  • Normality niet ok!

data_perpendicular<-Flights2 %>%  filter(Wind!="perpendicular")

model_perpendicular<-lmer(ForagingSpeed ~ WindSpeed_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=data_perpendicular)
summary(model_perpendicular)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ WindSpeed_KMI + (1 | NestID) + (1 | NestID:BaitID) +  
##     (1 | Individualcode)
##    Data: data_perpendicular
## 
## REML criterion at convergence: 259.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4499 -0.3729  0.0493  0.4221  3.3570 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.01943  0.1394  
##  NestID:BaitID  (Intercept) 1.12489  1.0606  
##  NestID         (Intercept) 0.86845  0.9319  
##  Residual                   0.28222  0.5312  
## Number of obs: 98, groups:  Individualcode, 47; NestID:BaitID, 40; NestID, 27
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)    3.41766    0.36128 50.36655   9.460 9.37e-13 ***
## WindSpeed_KMI -0.09641    0.06997 48.59291  -1.378    0.175    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.673
anova(model_perpendicular, ddf="Satterthwaite", type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## WindSpeed_KMI 0.53583 0.53583     1 48.593  1.8986 0.1745
res<-residuals(model_perpendicular)
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.95086, p-value = 0.001085
qqnorm(res)
qqline(res)

ForagingSpeed vs WindSpeed voor Upwind dataset

  • Hoogsignificant

  • Normality ok!

data_upwind<-Flights2 %>%  filter(Wind!="upwind")

model_upwind<-lmer(ForagingSpeed ~ WindSpeed_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=data_upwind)
summary(model_upwind)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ WindSpeed_KMI + (1 | NestID) + (1 | NestID:BaitID) +  
##     (1 | Individualcode)
##    Data: data_upwind
## 
## REML criterion at convergence: 432.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6138 -0.4276  0.1099  0.4081  2.2258 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1331   0.3648  
##  NestID:BaitID  (Intercept) 0.8066   0.8981  
##  NestID         (Intercept) 0.5930   0.7700  
##  Residual                   0.3105   0.5572  
## Number of obs: 174, groups:  Individualcode, 71; NestID:BaitID, 56; NestID, 28
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     4.13221    0.27666  64.34513  14.936  < 2e-16 ***
## WindSpeed_KMI  -0.22679    0.04392 122.09885  -5.163 9.55e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.631
anova(model_upwind, ddf="Satterthwaite", type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## WindSpeed_KMI 8.2778  8.2778     1 122.1  26.659 9.555e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res<-residuals(model_upwind)
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.98133, p-value = 0.01951
qqnorm(res)
qqline(res)

Model Selection

fullmodel<-lmer(ForagingSpeed ~ Traject100m + WindSpeed_KMI + Temperature_KMI + Cloudcoverage_KMI + Weight_ind + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
## boundary (singular) fit: see help('isSingular')
summary(fullmodel)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ Traject100m + WindSpeed_KMI + Temperature_KMI +  
##     Cloudcoverage_KMI + Weight_ind + (1 | NestID) + (1 | NestID:BaitID) +  
##     (1 | Individualcode)
##    Data: Flights2
## 
## REML criterion at convergence: 132.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4288 -0.5430  0.0748  0.5100  1.6174 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.07951  0.282   
##  NestID:BaitID  (Intercept) 0.66749  0.817   
##  NestID         (Intercept) 0.00000  0.000   
##  Residual                   0.25708  0.507   
## Number of obs: 71, groups:  Individualcode, 16; NestID:BaitID, 12; NestID, 7
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)        5.24098    2.24293 19.02657   2.337   0.0305 *
## Traject100m       -5.85766    2.63126  8.90757  -2.226   0.0533 .
## WindSpeed_KMI     -0.15017    0.20012 15.15305  -0.750   0.4645  
## Temperature_KMI   -0.05689    0.09633 21.97521  -0.591   0.5608  
## Cloudcoverage_KMI  0.59392    0.44256 38.57798   1.342   0.1874  
## Weight_ind        -1.12665    3.95322 11.66851  -0.285   0.7806  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 WS_KMI Tm_KMI Cl_KMI
## Traject100m -0.495                            
## WindSpd_KMI  0.684 -0.366                     
## Temprtr_KMI -0.809  0.309 -0.857              
## Cldcvrg_KMI -0.366  0.278 -0.308  0.397       
## Weight_ind  -0.794  0.146 -0.403  0.418  0.034
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
model1<-lmer(ForagingSpeed ~ Traject100m + WindSpeed_KMI + Temperature_KMI + Cloudcoverage_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ Traject100m + WindSpeed_KMI + Temperature_KMI +  
##     Cloudcoverage_KMI + (1 | NestID) + (1 | NestID:BaitID) +  
##     (1 | Individualcode)
##    Data: Flights2
## 
## REML criterion at convergence: 469.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.44419 -0.54207  0.03008  0.46341  2.91374 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1602   0.4002  
##  NestID:BaitID  (Intercept) 0.9180   0.9581  
##  NestID         (Intercept) 0.2815   0.5306  
##  Residual                   0.3234   0.5687  
## Number of obs: 191, groups:  Individualcode, 74; NestID:BaitID, 52; NestID, 23
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         4.06254    0.74726  64.14701   5.437 9.02e-07 ***
## Traject100m        -4.71721    1.54679  28.71297  -3.050 0.004887 ** 
## WindSpeed_KMI      -0.20065    0.05343 143.30993  -3.755 0.000251 ***
## Temperature_KMI     0.02873    0.03580  78.77035   0.803 0.424623    
## Cloudcoverage_KMI  -0.17268    0.21176 142.63991  -0.815 0.416163    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 WS_KMI Tm_KMI
## Traject100m -0.632                     
## WindSpd_KMI  0.385 -0.266              
## Temprtr_KMI -0.897  0.428 -0.593       
## Cldcvrg_KMI -0.343  0.097 -0.336  0.322
model2<-lmer(ForagingSpeed ~ Traject100m + WindSpeed_KMI + Cloudcoverage_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ Traject100m + WindSpeed_KMI + Cloudcoverage_KMI +  
##     (1 | NestID) + (1 | NestID:BaitID) + (1 | Individualcode)
##    Data: Flights2
## 
## REML criterion at convergence: 465.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.47342 -0.52496  0.04095  0.44893  2.92491 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1711   0.4136  
##  NestID:BaitID  (Intercept) 0.9100   0.9539  
##  NestID         (Intercept) 0.2484   0.4984  
##  Residual                   0.3224   0.5678  
## Number of obs: 191, groups:  Individualcode, 74; NestID:BaitID, 52; NestID, 23
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         4.59898    0.32492  32.71481  14.154 1.67e-15 ***
## Traject100m        -5.27920    1.37511  28.50382  -3.839 0.000632 ***
## WindSpeed_KMI      -0.17545    0.04319 132.44886  -4.063 8.27e-05 ***
## Cloudcoverage_KMI  -0.22488    0.20110 140.07082  -1.118 0.265372    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 WS_KMI
## Traject100m -0.620              
## WindSpd_KMI -0.422 -0.017       
## Cldcvrg_KMI -0.131 -0.047 -0.191
model3<-lmer(ForagingSpeed ~ Traject100m + WindSpeed_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ Traject100m + WindSpeed_KMI + (1 | NestID) +  
##     (1 | NestID:BaitID) + (1 | Individualcode)
##    Data: Flights2
## 
## REML criterion at convergence: 557.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.56525 -0.54882  0.06406  0.46203  3.02412 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1638   0.4047  
##  NestID:BaitID  (Intercept) 0.7251   0.8515  
##  NestID         (Intercept) 0.3533   0.5944  
##  Residual                   0.3176   0.5636  
## Number of obs: 230, groups:  Individualcode, 91; NestID:BaitID, 68; NestID, 34
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)     4.6054     0.2843  63.4475  16.197  < 2e-16 ***
## Traject100m    -5.1279     1.2857  41.4125  -3.989 0.000265 ***
## WindSpeed_KMI  -0.1895     0.0394 169.3697  -4.810 3.31e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100
## Traject100m -0.606       
## WindSpd_KMI -0.517 -0.011
library(cAIC4)
## Loading required package: stats4
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
## 
##     lmList
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:plotly':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
library(AICcmodavg)
## 
## Attaching package: 'AICcmodavg'
## The following object is masked from 'package:lme4':
## 
##     checkConv
AIC(fullmodel)
## [1] 152.8988
AIC(model1)
## [1] 487.598
AIC(model2)
## [1] 481.3893
AIC(model3)
## [1] 571.1946